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binaries in PTA data. We model the signal as a collection of circular monochromatic binaries, each 
characterized by three free parameters: two angles defining the sky location, and the frequency. We 
marginalize over all other source parameters and we apply an efficient multi-search genetic algo- 
rithm to maximize the likelihood function and look for sources in synthetic datasets. On datasets 
characterized by white Gaussian noise plus few injected sources with signal-to-noise ratio (SNR) 
in the range 10-60, our search algorithm performs well, recovering all the injections with no false 
positives. Individual source SNRs are estimated within few % of the injected values, sky locations 
are recovered within few degrees, and frequencies are determined with sub-Fourier bin precision. 
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I. INTRODUCTION 

Precision timing of millisecond pulsars provides a 
unique opportunity to get the very first low-frequency 
gravitational wave (GW) detection. This prospect is at- 
tracting the attention of the wider astrophysical commu- 
nity, causing a recent boost of activity in the field. The 
Parkes Pulsar Timing Array (PPTA, IH), the European 
Pulsar Timing Array (EPTA, [|) and the North Amer- 
ican Nanohertz Observatory for Gravitational Waves 
(NANOGrav, Q), joining together in the International 
Pulsar Timing Array (IPTA, are collecting data 

and improving their sensitivity in the frequency range 
of ~ 10~^ — 10~^ Hz. In the coming years, the Chinese 
five hundred meter aperture spherical telescope (FAST, 
0) and the planned Square Kilometer Array (SKA, Q) 
will provide a major leap in sensitivity. Current sur- 
veys arc already placing interesting upper limit on the 
level of a putative GW background [1, Q , skimming the 
range predicted by state of the art models of MBH evo- 
lution [1^ |Tl| . Within the next few years, the combined 
IPTA data might either result in a first detection, or start 
placing interesting limits on the MBH binary formation 
efficiency in massive galaxies. 

The detection principle is very simple: GWs affect the 
propagation of radio signals from the pulsar to the re- 
ceiver on Earth, leaving a characteristic fingerprint in 
the time of arrival of the radio pulses (e.g., [H, [l3|). 
Such fingerprint depends on the properties of the under- 
lying cosmological population of inspiralling binaries, and 
will consists of a superposition of quasi-monochromatic 
waves, similar to the white dwarf- white dwarf foreground 
(e.g., in the mHz window relevant to space based 



interfcrometry 15|. This signal has generally been re- 
garded as a stochastic background, and data analysis 
techniques has been developed accordingly [1, H, ■ 
The actual expected signal, however, is far from being 
isotropically distributed in the sky, with just few sources 
dominating the power at each frequency [13, [HI H^l ■ The 
possibility of resolving an individual source offer appeal- 
ing astrophysical prospects, and PTA capabilities on this 
front were also investigated in details by many authors 

What is still missing, is a detailed study of what kind 
of information a PTA can extract out of a complex su- 
perposition of multiple sources. Is the signal going to be 
similar to a confusion noise? Can we resolve individual 
sources? How many of them? Can we locate them in the 
sky, and to what level of accuracy? All these questions 
are of great interest for the astrophysical community; pre- 
cise sky localization of individual sources will allow the 
efficient search for electromagnetic counterparts [l^, [111 , 
opening the new horizon of multimessenger astronomy. 



This is a second in a series of paper devoted to the 
exploration of the PTA potential of resolving multiple 
GW sources. In [l| (hereinafter Paper I), wc demon- 
strated PTA efficiency in disentangling monochromatic 
sources at the same frequency. The key idea is to esti- 
mate the likelihood that a certain number of sources with 
certain parameters are present in the data. We developed 
a formalism that allowed us to maximize analytically the 
likelihood over the extrinsic source parameters, restrict- 
ing the search to the source sky location only (2x7V pa- 
rameters, where is the number of GW sources in the 
template). There, we did not implement any proper algo- 
rithm to search over the parameters space, and we made 
a lot of simplifying assumptions, suitable to a first, ex- 
ploratory investigation. 

Our aim is to implement a proper search algorithm, 
progressively relaxing our limiting assumptions to de- 
velop a detection pipeline able to handle the whole com- 
plexity of a realistic datasct. We start in this paper with 
two major steps: (i) we extend our mathematical for- 
malism to include frequency scan and (ii) we present an 
upg raded version of the genetic algorithm employed by 
j29j | in the LISA mock data Challenge [10, [lH specifi- 
cally developed to search for a global maximum on the 
multimodal likelihood surface embedded in the multidi- 
mensional parameter space. We have found (similarly 
to the mock LISA data challenge) that the genetic algo- 
rithm (GA) is very efficient in finding the correct number 
of sources and their parameters. 

The paper is organized as follows. In Section 2, we spell 
out our main assumptions and in Section 3 we present 
the genetic algorithm and its feature. The datasets used 
to test the algorithm are detailed in Section 4, and the 
algorithm performances and results arc presented in Sec- 
tion 5. In Section 6 we draw our conclusion and discuss 
improvements we will present in future work. 



II. DETECTION STRATEGY, EXTENSION TO 
FREQUENCY SEARCH 

The main purpose of this paper is to extend our for- 
malism to include search in frequency and to implement 
a proper search algorithm to identify maxima in the like- 
lihood. Accordingly, we relax number 1, 3 and 8 from 
the limitations and assumptions described in Section II 
of Paper I, i.e.: 

1. wc consider only datasets with noise; 

2. wc inject sources at different frequencies; 

3. we implement a proper search algorithm to maxi- 
mize the likelihood. 



A. Choice of the template 

In computing the likelihood function, we consider 
monochromatic GW sources assuming that the orbital 
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frequency does not chan ge a ppreciably over the obser- 
vation period (see, e.g., [22|) which we took to be 10 
years. Each GW signal is therefore characterized by 
seven parameters only: the overall signal amplitude {-4}, 
the source frequency and phase {/, $o}j and the angles 
defining its location in the sky {(j),9}, inclination {(,}, and 
polarization {ip}- Contrary to Paper I, we do not fix all 
the systems at the same frequency but wc consider the 
unknown frequencies of the sources as additional search 
parameters. 

Wc do not use the full response of the pulsar-Earth 
detector, which is combined out of the "Earth term" and 
the "Pulsar term" (see, e.g., [22j) (from now on we drop 
the quotes and use those notions as a jargon), but we con- 
struct a signal template to match the Earth term only. 
There arc several reasons to 'drop' the pulsar term in the 
analysis. In all pulsars, the Earth terms add-up coher- 
ently: they all have the same frequency and phase, and 
the amplitude of the signal in the residuals depends on 
the relative position of the pulsar and the GW source 
on the sky. Conversely, the pulsar terms are in general 
incoherent: they usually appear at different frequencies, 
and the phase and amplitude of the signal depends not 
only on the position of the source relative to the pulsar, 
but also on the distance to the pulsar which is usually 
poorly known (in most of the cases to ^ 10% precision). 
Even if we assume we know the pulsar distance exactly, 
the pulsar term carries the imprint of the binary system 
emitting as it was a time At = L(l -|- k.n) in the past as 
compared to the Earth term; where L is the Earth-pulsar 
distance and k and n are the unit vectors pointing to the 
sky location of the source and the pulsar respectively. 
This means that to connect pulsar and Earth terms we 
need to know the evolution of the binary system for At 
which is typically 10^ — 10^ years. Even assuming pure 
GW evolution, the prediction of the signal at the pulsar 
term will be affected by spin-orbit coupling precession 
js^ . or non negligible eccentricity (which is very likely 
for broad binaries (sec, e.g., [33l.[34|). In addition to this, 
the assumption that the binary evolution is driven by 
GW back reaction only could not hold for such system. 
Widely separated MBH binaries could still dynamically 
interact with the surrounding gas and / or stellar environ- 
ment, which might significantly affect (and sometimes 
dominate) their orbital evolution In other words, 

the inclusion of the pulsar term is always model depen- 
dent, and we try to avoid it. Considering the Earth term 
only also has some drawbacks. If the systems evolve only 
under GW radiation reaction, then relatively light and 
wide binaries might not have evolved much (less than a 
Fourier bin) over the pulsar-Earth light travel time; in 
this case the pulsar and Earth terms appear effectively 
at the same frequency and we need to take that into ac- 
count. Moreover, pulsar terms from different GW sources 
will inevitably overlap in frequency with the Earth term 
of a given source, creating a spurious contribution to the 
signal. Here wc inject in the data only the Earth portion 



of the signal, and we delegate problems related to pulsar 
terms to future work. 



B. Likelihood function and detection statistics 

The details of the detection statistics were outlined in 
Paper I, we briefly summarize here the main points and 
describe the extension of the formalism to sources with 
different frequencies. As justified in the previous sec- 
tion, we use the matched filtering technique assuming the 
Earth term as a template. The mathematical description 
of the signal template for an individual source as a func- 
tion of the parameters A = {A, f,^o, tjip} is given 
by equations (11)-(16) of Paper I. The log-likelihood ra- 
tio (likelihood that a dataset Xa (t) contains a GW signal 
ra(t] A) over the likelihood that it is pure noise), is 

logA =< Xclr^ > < r„|r„ >, (1) 

where the subscript a corresponds to a given pulsar and 
Tq, = is the expected Earth term in the data. Wc 
neglect here all possible stochastic GW signals and look 
for individual binaries standing above a putative unre- 
solved background only. The inner product appearing in 
equation ((!]) is defined as 

TV 

<^k>=^^E^(^'MiO, (2) 

where N is the number of points in the time series, Tq is 
the observation time, and, >S'(/) is one-sided noise power 
spectral density which we assume to be white Gaussian. 
Equation ([2]) is the discrete version of the inner product 
used in Paper I; it has the advantage to be applicable to 
unevenly sampled data, which will be the case in reality. 
It was shown in Paper I, that the GW signal imprinted 
in the data by each individual source can be written as 

4 

= EaO)/»0), (3) 

where are time dependent functions that include the 
parameters we want to search for, while ajj) are constants 
over the observation period. Expressions for hf^j^ and a^j) 
are the same as given in Paper I. Here, we search over 
sky location and frequency / of each GW signal, but we 
assume / to be constant over the time of observation, 
therefore we can keep it either as part of a(j) (as in the 
Paper I) or we can include it in h(^j^ . In practice we 
decompose the signal in hi^j){9,4>, f) and oq) {A, $0 , t, ■(/') 
(therefore shifting the / dependence from oq) to 
compared to the expressions given in Paper I). 

Wc can then maximize the likelihood ratio over the a^j) 
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constants for each GW source analytically: 



glog(A) 
da 



= 0, ^a(fc)=A/fc/Xj, 



kj 
--] 

jk ^^j' 



(4) 
(5) 



where 



The statistical properties of J>, arc investigated in details 
in [2^. In presence of Ng GW sources in the template, 
the coefficients a^^) are represented by a 4 x A^^ array, Xj 
is also a 4 X array, while the M-matrix is a ANg x 
4iVs, 2-D matrix. The matrix can be decomposed in Ng 
4x4 row- and column-matrices, each corresponding to 
the cross terms between the I-th and J-th GW sources: 



M 



I, J 



E 



Q"iig^ v"xli Q"xli v"xli 

U'-'Xli Q'-'Xli U'-'Xli Q^'^X'J 

\ Q^^xii y"^ Q"^^ y"^li ) 



(7) 



where 



V 



IJ 



(8) 



and cosine combinations are given by: 
To Jo 

= sinc(A0) — sinc(I]0) 



rlJ 



cos{uj^t) cos{u}'^t) dt 



To Jo 
sinc(A(/)) -I- sinc(I]0) 



To 



sin(w^t) cos{uj"'t) dt 



(10) 

(11) 
(12) 
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— I CQs{uj t)sin{(jj t) dt 



(13) 



o Jo 

sm 1 — 1 smc ( — 1 - sm ( — 1 smc ( — 



where A0 = (cj-^ - w"')To and T,(p = {uj' + uj-')To. Note 
that the M matrices reduce to the expression given in 
equation (25) of Paper I when o;^ = w'. The I ^ J terms 
give beatings between two signals at different frequencies 
and they are usually smaller than the terms in the I = J 
matrices. We found that one can consider the sources ap- 
proximately at the same frequency if \fi—fj\ ~ (2/3)AF, 
where AF = I/Tq is the size of the Fourier frequency bin. 

We use jFe as detection statistic. Note that we can 
also estimate the relative contribution of each source as 
J"/ = \X({M-YY . Following we can express 
the relation between the expectation of the analytically 
maximized likelihood J> and the SNR as follow: 



'ms + SNR^ 



(14) 



To search for an individual source we use the same math- 
ematical framework assuming Ng = 1 , we refer to [26j for 
more details on the statistical properties of J-e in this 
latter case. 



and F^^ represent the decomposition of the antenna pat- 
tern given by equation (16) of Paper I. The X terms come 
from the inner products of the time dependent parts of 
/i^-j which, for each source / are cos- and sin-functions 
of phase (/)/ = lufji = ujit. We can evaluate those inner 
products analytically by using the integral representation 
adopted in Paper I; for example: 



III. MULTI-SEARCH GENETIC ALGORITHM: 
DESCRIPTION AND IMPLEMENTATION 

We search for the maximum of J^e with a modified 
version of the genetic algorithm (GA) described in poj . 
performing multiple searches in parallel. 



< ^ImKlD iFcYiF^Y < sin(0/)| sin(0^) > 
^{F^^F^Y^ I °sin(0z)sin(0^)dt 

-to Jo 

^ {F:Y{Fn'Tii ^ u"xig'. (9) 

The explicit form of the X integrals for all possible sine 



A. Genetic algorithm 

The GA is a method to perform global searches on 
large parameter spaces (optimization method) based 
on the natural selection principle. In nature, or- 
ganisms adapt themselves to their environment: the 
smartest/strongest/healthiest organisms are more likely 
to survive and participate in the breeding to produce the 
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Genetic algorithm 



GW search 



organism 
gene (of an organism) 
allele (of a gene) 
quality Q 
colony of organisms 
71-th generation 
(selection + breeding) + mutation 



template : signal from A'^s GW sources 
parameter (of a template) : 3 x A''^, 
bits (of the value of the parameter) 
Maximized Likelihood, i.e. F-statistic 
evolving group of templates 
the state of colony at n-th step of evolution 
w parameter space exploration strategy 



TABLE L Correspondence between GA and GW data analysis notions. 



offspring. These two processes, selection and breeding, 
are used in GAs to produce subsequent generations of 
organisms. Since the best organisms are more likely to 
participate in breeding, the new generation should be 
better (in which sense we will specify later) than the pre- 
vious one, or at least no worse. This procedure leads to 
an evolution of the organism population, just like in na- 
ture: the good qualities of the parents can be transferred 
to their offspring. In the biological world, besides these 
two basic operations, among every generation, there are 
always few individuals which have better characteristics 
to adapt to the environment, produced as a result of pos- 
itive mutations. By introducing new genotypes into the 
population, mutations can potentially improve the forth- 
coming generations and consequently accelerate the evo- 
lution towards a perfect adaptation to the environment. 

We apply these principles to the search for individual 
GWs in PTA data using equivalences described in tablcUl 
a template described by a set of parameters {6'/,0/, //} 
is one organism described by a set of genes; the of 
the template is the quality of the organism; the binary 
representation of a parameter by a set of bits is the rep- 
resentation of a gene by a set of allele. 

We start with a group of organisms (templates) chosen 
randomly (initial search) or constructed from the results 
of previous searches. We evaluate the quality of each or- 
ganism (J-"e). We select set of pairs (parents) based on 
their quality: organisms with better quality (templates 
with higher T^) are chosen more often than weak or- 
ganisms. We combine the genotypes of two parents to 
produce a child (we combine parameters of two chosen 
templates to produce a new one) . We impose the number 
of produced children to be equal to the number of parents 
(i.e., we keep the number of evolving organisms constant 
at each generation). Next, we allow with a certain prob- 
ability a random mutation in the children's genes (with 
some probability we randomly change the parameters of 
the new templates, exploring a larger area of the param- 
eter space). The parents are discarded and the resulting 
children form a new generation. We repeat the procedure 
until we reach a steady state (maximum in the quality), 
or up to a maximum number of generations. We keep 
only one generation active (one group of templates). 

We now turn to describe the tree steps used for making 



a new generation (replacing parents by children), speci- 
fying the used possibilities of tuning for each process. 

In the selection process, two parent organisms are se- 
lected. The probability to chose an organism is directly 
related to its quality J-g: the higher is Te of an organism, 
the higher is the probability to be selected. The relation 
between and the selection probability is moderated by 
a parameter called "temperature" . The higher is the tem- 
perature, the smaller is the difference in selection proba- 
bility among the organisms (for infinite temperature, all 
probabilities are equal); the lower is the temperature, the 
higher is the probability to select only the best organisms. 
The value of the temperature evolves during the search 
(similar to simulated annealing [1^): it starts at high 
temperature and then decreases and alternates between 
hot and cold phases: i.e., we allow at the beginning a 
large exploration range, by accepting good and bad or- 
ganisms, but later on we search only around the better 
ones, allowing some jumps (by alternating hot and cold 
phases). More details about the selection criteria can be 
found in sections III.C and IV.A.2 of [H]. 

In the breeding process, the parameters of the two se- 
lected templates are mixed by combining the bits of their 
binary representation, i.e. by combining the allele of the 
genes. More details on the breeding and code of the genes 
can be found in sections III.B and III.D of [l^. At each 
generation, we change the binary representation alter- 
nating between "standard binary" and "gray code". The 
breeding method used here is the " crossover one random 
point" which consists of taking the first bits from one 
parent and the last ones from the other ; the cross point 
is chosen randomly. 

In the mutation process, some of the bits are changed. 
The mutation rate is managed by a parameter called 
PMR (Probability Mutation Rate). At the beginning of 
a search, a gene is mutated with a probability described 
by the PMR. Then for mutating a gene, 8 bits of this 
gene are (over 20) randomly chosen and changed. After 
a certain number of generations (typically 300), the type 
of mutation is changed : each bit of each gene is directly 
mutated with a probability at PMR. The value of the 
PMR decreases during the search, starting around 0.5 - 
0.1 and ending around 0.1 - 0.01. Using these two types 
of mutation and the PMR evolution corresponds again 
to start with a large exploration, slowly shrinking to the 
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area around the best templates only at late generations. 
More details on the mutation are provided in sections 
III.E and IV.A.3 of 

In addition, we always reproduce the best organism 
between two generations (elitism - cloning of the best). 
This means that the algorithm always converges toward 
the best solution. Finally, we also use the local mutation 
described in section IV. B. 2 of (29l |. 

We typically use 50 organisms per generation and 1000 
generations. The run of one GA takes few minutes on a 
standard laptop (one Intel core at 2GIIz). Since the size 
of the parameters space increases with the number of 
sources in the template, the convergence speed decreases 
accordingly. The algorithm usually converges around the 
true solution in less than 400 generation for the highest 
SNR sources. One of the most interesting features of the 
GA is its efficiency in finding maxima in the J-e surface 
first (during the large exploration phase), and then in 
exploring them deeply to extract the global one (during 
the local exploration phase) . One G A run is usually suf- 
ficient to find most of the sources, but sometime it gets 
stuck on some local maximum. To overcome this prob- 
lem, wc run several GAs in parallel, as described in the 
next section. 



B. Multiple searches (MultiSearch) 

The GA described in the previous section provide 
the basis for a more general method called "multiple 
searches" (MS) algorithm. This method consists in run- 
ning several GAs in parallel with different properties and 
initial parameters. 

Wc take an initial population with parameters chosen 
randomly. Wc start a GA on this population, tuning 
the parameters to perform a large exploration. In the 
resulting population, we select only the best organism 
which are well separated. This means that the selected 
organisms have SNR > 97% SNReost and the distance in 
parameter space between two organisms is higher than 
a certain threshold chosen empirically after a number of 
tests: I cos{9j^i) — cos{9jj)\ > Ace = 0.1 , — (/i/^ l > 
= 20° and \fj,i — fj,j \ > A/ = 0.5nHz, where / refers 
to the source and i and j to the solutions. The selected 
solutions are called " modes" and this selection process is 
called " mode separation" . 

The next step is to start one GA on each mode, tuned 
for local exploration. The goal is to explore the vicinity 
of the mode to find the local highest value of Jg. The 
organisms of each GA, are allowed to explore only their 
mode neighborhood and are forbidden to go on the area 
of interest of other modes. The area of interest of a mode 
{cos 9i, (f>i,fi} is defined within [cos^i — Acg, cos 9i + Ace], 
[04 - A^,(t)i + A^] and [fi ~ Af, fi + Af]. In parallel to 
these 'mode GAs', we start another GA tuned for large 
exploration. We forbid the organisms of this GA to go 



on the areas of the modes. The aim of this GA is to find 
new modes (overlooked in previous searches) , if there are 
any left (it can also give a null result). 

At the end of this step, all the solutions are grouped 
together and we apply the " mode separation" to identify 
"modes". Then we iterate the procedure by restarting 
several 'local' GAs. 

In the long run, this method, as other stochastic meth- 
ods (e.g., Markov Chain Monte Carlo methods), is guar- 
anteed to converge to the global maximum. However, 
there is no way to exactly know a priori how fast it will 
do so, and one has to decide when to stop it, being some- 
how confident that the best solution has been found. We 
usually do 2 to 5 iterations of the procedure outlined 
above before stopping. The number of modes N^o^cs 
found increases with the number of iterations. Since we 
are running A^modcs + 1 GAs at each iteration, the first 
one takes just few minutes (one initial exploratory GA 
run), the second one can take up to an hour (depend- 
ing on the number of modes found in the first iteration), 
and the later ones up to few hours. In total, we run be- 
tween 50 to 300 GAs for a search. The correct solution 
is usually found after 2 iterations (i.e. about one hour). 
As a pseudo-test for convergence, we run several times 
(typically 10) our MS-GA code, with different initial con- 
ditions. If all the run give almost the same results, we 
claim convergence. 



IV. DESCRIPTION OF THE TEST DATASETS 

The genetic algorithm described in the previous sec- 
tion was used to analyse four blind datasets, which we 
describe here in detail. Each dataset consists of a col- 
lection of time series representing the residuals obtained 
by timing an ensemble of millisecond pulsars (MSPs). 
In all datasets, MSPs are placed randomly in the celes- 
tial sphere, each time series consists of 523 equally sam- 
pled datapoints over a total observing time of 10 years 
(one datapoint every two weeks), and the noise is as- 
sumed to be white Gaussian. The injected sources were 
all equal mass, circular, non spinning binaries with chirp 
mass of 10^ Mq, placed at the same rcdshift (distance), 
but with sky location, inclination, polarization and ini- 
tial phase drawn randomly, resulting in a range of signal 
strengths. The redshift was chosen to produce the de- 
sired SNR range in the datasets (see below). The fre- 
quency was drawn from a random distribution in the 
range 10~* — 10~^Hz. Sources were evolved according to 
equation of motion accurate to 3.5 Post Newtonian or- 
der in phase evolution 13711 . and gravitational waveforms 
were generated following j22l| (see also [3l,|3§|)- The final 
residual injected in the datasets were obtained by time 
integration of the waveforms (see equations (8) and (9) of 
Paper I). Note that the injected data are quite different 
from the adopted circular, non evolving monochromatic 
templates we use in the search; we are therefore mimick- 
ing the (likely) situation in which the template does not 
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FIG. 1. Sample of simulated timeseries. In each panel, the red dashed curve is the injected signal, where the blue jagged line 
represent the total raw dataset including signal plus white Gaussian noise. Left panel: pulsars extracted from DatasetS; sources 
are injected in white Gaussian noise with rms a — 50ns. Right panel: pulsars extracted from Dataset4; here each pulsar has a 
different noise level, as labelled in each panel. 



perfectly match the signal. Especially at high frequency, 
there might be a non-negligible evolution of the source 
frequency over 10 years, possibly introducing a bias in 
our source recovery. We will quantify this effect in our 
results. As in Paper I, we considered the Earth term only 
(issues related to pulsar terms will be explored in the next 
paper). Here follow the details of the four datasets: 

• Datasetl: 30 MSP; rms noise 50ns in each pulsar; 
5 binaries at z — 0.01, with individual SNR in the 
range ^ 30 — 60; 

• Dataset2: 30 MSP; rms noise 50ns in each pulsar; 
4 binaries at z = 0.02, with individual SNR in the 
range ~ 15 — 55; 

• DatasetS: 50 MSP; rms noise 50ns in each pulsar; 
8 binaries at z — 0.03, with individual SNR in the 
range 10 — 40; 

• Dataset4: 50 MSP; rms noise of each pulsar ran- 
domly drawn in the range 30 — 200ns; 3 binaries 
at z = 0.01, with individual SNR in the range 
- 30 - 40. 

Datasets are in order of increasing complexity (more 
sources, lower SNR). In the last dataset we tested the 



algorithm performance when combining time series with 
different noise levels. Sample time series extracted from 
Dataset3 and Dataset4 are visualized in figure [TJ where 
we can appreciate the variety of imprints depending on 
the pulsar location in the sky relative to each individual 
source. 



V. RESULTS AND DISCUSSION 

The datasets were generated separately by A. Sesana 
and were blindly analyzed by A. Petiteau and S. Babak. 
The MS-GA was applied to all datasets, adding sources 
one by one to the template. By doing this, we could 
test the effectiveness of the code in determining both the 
number of sources in the dataset and their sky location. 
A summary of the results is given in table |lll 

For each dataset we evolved several colonies of or- 
ganisms assuming Ns = 1,2,3,.... in the template, we 
computed the SNR of the best organism at the end of 
each search, and track its evolution with Ng. Results 
are shown in the left panel of figure [2] for Dataset3. The 
maximum SNR steadily increases by adding sources up to 
Ns = 8. Adding a ninth source to the template does not 
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FIG. 2. Performance of the MS-GA in finding the number of sources. Left paneh signal SNR as a function of the number of 
sources Ns assumed in the template. Right panel: source localization in the sky (sky map) and in the frequency-SNR space 
(rightmost panel). Blue crosses (x) correspond to injected values and black dots to the position of the MSPs forming the array. 
Red marks represent sources l-to-8 found by all the organisms with SNR^^t > 99% SNR^osti while green circles represent the 
9th source. Note that this latter one does not have a defined position, and typically has SNR< 5. 



TABLE IL Recovered and (injected) parameter values of all the simulated sources in each dataset. The last column represent 
the sky offset of the recovered sources with respect to the injection (see text for details). 





SNR 


.f[ns] 


6l[rad] 


4>[iaxi] 


Ae[dcg] 




60.70 


(61.11) 


56.6 (56.4) 


1.249 


(1.237) 


2.604 


(2.601) 


0.706 




45.85 


(42.48) 


38.0 (38.0) 


1.750 


(1.748) 


3.765 


(3.764) 


0.127 


Datasetl 


43.71 


(40.43) 


36.4 (36.4) 


1.555 


(1.529) 


1.722 


(1.712) 


1.596 




35.67 


(36.14) 


53.7 (53.5) 


0.537 


(0.534) 


5.522 


(5.451) 


2.085 




32.27 


(31.33) 


48.3 (48.0) 


1.286 


(1.295) 


5.144 


(5.123) 


1.266 




54.64 


(54.07) 


18.88 (18.9) 


1.774 


(1.774) 


3.839 


(3.841) 


0.112 


Dataset2 


48.01 


(47.24) 


11.25 (11.3) 


1.870 


(1.858) 


5.720 


(5.718) 


0.696 




13.64 


(13.05) 


77.42 (76.5) 


0.617 


(0.651) 


6.158 


(6.115) 


2.434 




12.23 


(12.78) 


57.19 (57.0) 


1.613 


(1.549) 


6.050 


(6.048) 


3.669 




44.91 


(42.99) 


19.33 (19.3) 


0.474 


(0.468) 


1.450 


(1.454) 


0.359 




37.39 


(37.72) 


25.42 (25.4) 


0.883 


(0.878) 


2.733 


(2.749) 


0.763 




26.02 


(27.09) 


13.21 (13.2) 


1.769 


(1.764) 


5.078 


(5.087) 


0.581 


DatasetS 


20.19 


(20.88) 


83.42 (82.4) 


0.689 


(0.668) 


4.133 


(4.162) 


1.593 




19.67 


(18.51) 


39.79 (39.8) 


0.541 


(0.509) 


0.386 


(0.429) 


2.211 




17.27 


(16.59) 


33.16 (33.1) 


1.381 


(1.397) 


3.621 


(3.693) 


4.160 




13.07 


(13.19) 


73.83 (73.0) 


1.534 


(1.536) 


5.054 


(5.078) 


1.379 




10.66 


(11.51) 


82.75 (81.8) 


0.809 


(0.864) 


6.182 


(6.085) 


5.192 




42.73 


(43.92) 


98.2 (96.3) 


2.028 


(2.043) 


0.977 


(0.961) 


1.200 


Dataset4 


28.62 


(29.28) 


91.5 (90.1) 


2.655 


(2.661) 


1.174 


(1.121) 


1.454 




27.56 


(28.28) 


48.2 (48.1) 


1.231 


(1.245) 


5.774 


(5.769) 


0.827 



significantly improve the match with the data, indicat- 
ing that the dataset is best described by an eight-source 
model; in fact there were eight sources in DatasetS. The 
algorithm identified the correct number of sources in all 
datasets. We stress here that all the injected sources had 
SNR> 10; high enough to be dug out of the noise. In 
presence of many low SNR sources, we do not expect 
any search algorithm to recover the correct number of 



binaries, but only to identify the brightest ones. We will 
address this 'confusion problem' in a future paper. A 
complementary view of this result is given in the right 
panel of figure [2] We ran several GAs using nine-source 
colonics of organisms, we identified the solution (organ- 
ism) with highest SNR (SNRbost) and we stored all the 
organism having SNR^^^t > 99% SNR^^.^. The figure 
show the location in the sky of all sources found in all 
these 'best solutions'. The location of sources l-to-8 does 
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not change much for different solutions, and it is gener- 
ally consistent with the true (blue crosses) locations of 
the injected sources. Conversely, the 9th source (green 
circles) is extremely scattered around the sky. Moreover, 
the frequency-SNR plot at the extreme right shows that 
the individual SNR of those 9th sources are almost al- 
ways < 5, compatible with noise fluctuations. 

Having tested the code effectiveness in finding the 
number of sources present in the data, we turn now 
to the description of the results obtained on the indi- 
vidual datascts. Best solutions (those with SNR^q^ > 
99.5% SNR^^^J for Datasetl are shown in the top panel 
of figure [31 All the five injected sources were found at 
approximately the right sky location, with the right fre- 
quency and SNR. Our GA is designed to find the modes 
corresponding to maxima in the likelihood function, but 
not to explore the exact shape of the likelihood func- 
tion around those modes. The lack of parameter space 
exploration around the maxima prevents us to attach 
fully meaningful errors to our best solutions. We plan 
to include systematic exploration of the maxima in fu- 
ture work, here we just estimate the sky location error 
as the angular offset between the best solution and the 
injected signal. This is defined as AO — arccos(nt • fir), 
where fit and fir are the unit vectors defining the true 
sky location of the sources and the recovered value re- 
spectively. This is reported in the last column of table 
im All sources in Datasetl are offset by less than 2deg. 
Results for Dataset2 are shown in the central panel of 
figure [3l Again, we see that all sources are correctly 
identified, despite two of them having SNR just above 
10. Sky location offsets A8 are less than Ideg for the 
two low frequency bright sources, but degrade to '--^Sdeg 
for the high frequency, faint ones. DatasetS was the rich- 
est of all, with eight injected sources. Best solutions are 
shown in figure HI for different SNR threshold, to give a 
sense of 'how fast' points cluster toward the maximum 
of the likelihood. Also in this case, all sources are well 
located in the sky, with brighter sources located better. 
Looking at table |lTl we notice that we tend to oversti- 
mate the frequencies of sources above 60 nHz. This is 
because at such high /, the 10^ Mq chirp mass binaries 
injected in the data chirp significantly over the 10 year 
duration of the observations. Since we are matching the 
signal with non-evolving monochromatic templates, the 
estimated frequency is higher than the one injected at 
the beginning of the observation. However, such high 
frequency sources seem to have sky location offsets simi- 
lar (maybe a little worse) to their low frequency counter- 
parts with comparable SNR; i.e., the frequency mismatch 
between the signal and the template does not seem to 
corrupt the sky location performance of the search. The 
same effect is seen in Dataset4 (bottom panel of figure 
13] and table . Also in this case we find source offsets 
within ~ldeg of their true position, but we give a cou- 
ple of extra nHz to the high frequency sources. Different 
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FIG. 3. Best solutions for Datasetl (top panel), Dataset2 
(central panel) and Dataset4 (bottom panel). In each panel all 
solutions with SNR?ot > 99.5% SNR^ 

est are shown. Symbols 
have the same meaning as in the right panel of Figure (2] All 
recovered sources are color-coded according to the rightmost 
scale, based to the total SNR of the solution they belong to. 



noise levels in the pulsars do not affect the performance 
of our search algorithm. 

Overall, our MS-GA performed well on all datasets, 
recovering all the injected sources without returning any 
false positive. The parameters of the recovered sources 
well matched the injections with: (i) sky location offsets 



10 



src 1 
src 2 
src 3 
src 4 



src 5 
src 6 
src 7 



-1*' 


:■;«■■■■•, 

«^ '^^^ 













src 1 □ 
src 2 Q 
src 3 ^ 



src 5 
src 6 
src 7 



psr ■ 
true X 



src 4 V src 8 





© ■ - 


■, 











src 1 
src 2 
src 3 
src 4 



src 5 
src 6 
src 7 
src 8 



psr 
true 



> 








- 








































i . 




j 

i 


! 






















I 













20 40 60 80 100 
Frequency (nHz) 



-5 


1 














































1 - 









































20 40 60 80 100 
Frequency (nHz) 



> 


3 




























< 


















i - 




> 


< 














i 



74 

73.9 

73.8 

73.7 

73.6 

73.5 

73.4 

73.3 

73.2 

74 

73.9 

73.8 

73.7 

73.6 

73.5 

73.4 

73.3 

73.2 

74 

73.9 

73.8 

73.7 

73.6 

73.5 

73.4 

73.3 

73.2 



20 40 60 80 100 
Frequency (nHz) 



FIG. 4. Best solutions for DatasetS. The top panel shows 
all solutions with SNRjot > 99% SNR^^^t, the central panel 
all solutions with SNR^^t > 99.5% SNRbost and the bottom 
panel all solutions with SNR?ot > 99.8% SNRtcst- Symbols 
have the same meaning as in the right panel of Figure [J] All 
recovered sources are color-coded according to the rightmost 
scale, based to the total SNR of the solution they belong to. 



however take sky position offsets as a proxy of the sky 
localization accuracy. In fact, offsets shown in the last 
column of table HIl scale (with a large scatter) with the 
inverse of the SNR. This has to be expected: an off- 
set scaling with 1/SNR implies an area of uncertainty 
scaling with 1/SNR^, in agreement with theoretical ex- 
pectations. If we approximate the crrorbox in the sky 
as AO « 7r[A0]^, we get values in the range 10-70 deg^ 
for sources with SNR in the range 11-13. This is broadly 
consistent with [22| who estimated an average sky loca- 
tion accuracy of Af2 sa 50deg^ for a source observed by 
an array of 50 pulsars, randomly located in the sky, with 
total SNR= 10 (in the Earth term). Another interesting 
fact is the frequency mismatch for sources approaching 
10~^ Hz, caused by their frequency evolution over the ob- 
serving time (10 years). This means that, in principle, for 
such sources, we can measure the frequency drift /, i.e., 
the chirp rate. In fact, it takes only an extra parameter in 
the template, with some extra computational cost. The 
measure of / breaks the chirp mass/luminosity distance 
degeneracy in the source amplitude, allowing for a direct 
measurement of the source luminosity distance. This, 
ultimately, will narrow down significantly the number of 
candidate electromagnetic counterparts in the source sky 
errorbox, facilitating a positive identification. However, 
evolution on such short timescales is detectable only for 
very massive (A4 ^ 10^ Mq, as the systems injected in 
the data) binaries emitting at frequencies higher than 
'--^ 7 X 10~*Hz. Intrigued by this possibility, we checked 
how likely is to find such extreme systems in realistic 
populations of MBH binaries in the Universe. We took 
the models investigated by [Tl| and computed the av- 
erage number of expected sources with > 10^ M0 
and / > 7 X 10^^. Depending on the adopted MBH 
mass-bulge relation and on the accretion implementation 
(see [m for details), we found average number of sources 
ranging from 10"'^ to 0.04, i.e., there is less than 5% 
chance to have a such bright high frequency source in the 
sky. If we relax the mass requirement to > 5x10® Mq , 
figures grow to 0.01-0.4. To properly quantify the prob- 
ability of measuring /, one should estimate its minimum 
measurable value for a given array, and then select in 
the MBH binary population all the sources occupying 
the portion of the chirp mass-frequency parameter space 
compatible with such value. The crude figures estimated 
here indicate that / measurements using the Earth term 
only should be unlikely. 



VI. CONCLUSIONS 



less than few degrees, (ii) individual source SNR estima- 
tions within few % of the true ones, and (iii) sub- Fourier 
bin frequency accuracy (sometimes within 0.1 nHz for 
low frequency sources). Without a complete exploration 
of the likelihood function around the maxima, it is dif- 
ficult to asses proper errors on the parameters. We can 



This is a second in a series of paper devoted to the 
exploration of the PTA potential of resolving multiple 
GW sources. In Paper I we addressed basic issues like the 
number of sources per frequency bin that can be resolved 
by an array of N pulsars, demonstrating our findings with 
primitive searches on several (mostly noiseless) synthetic 
datasets. Here we pushed our analysis a bit further by (i) 
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extending the mathematical formulation of the likelihood 
function to include the source frequencies as additional 
free parameters and by (ii) implementing a multi-search 
genetic algorithm to efficiently find the maximum of the 
likelihood function. 

Wc constructed synthetic datascts consisting of collec- 
tions of time series representing the residuals obtained 
by timing an ensemble of MSPs. MSPs were placed ran- 
domly in the sky, each time series consisted of 523 equally 
sampled datapoints over an observing time of 10 years 
(one datapoint every two weeks), and the noise in the 
data was assumed to be white Gaussian. In each dataset 
we injected an unknown number Ns of sources with ran- 
dom parameters and individual SNR > 10 and we apply 
our multi-search genetic algorithm to search for their sky 
location and frequency. Note that we assumed circular 
monochromatic sources in our template, but we allowed 
for full PN evolution of the injected sources. By doing 
so, we placed ourselves in the (likely) situation in which 
the theoretical model of the signal does not perfectly rep- 
resent its real nature, and we explored the consequences 
of this mismatch. 

Our main results can be summarized as follows: 

• the MS-GA generally converged to the true maxi- 
mum of the likelihood function in 2-to-5 iterations 
(few hours on one core at 2GHz); 

• the MS-GA successfully identified all the injected 
sources in all datascts. No false positive were found; 

• the search on all source parameters was success- 
ful: inferred sky locations were offset by less than 
few degrees, individual source SNR estimations 
matched the injections within few %, and frequen- 
cies were determined with sub- Fourier bin precision 
(most of the times to better than O.lnHz); 

• sky location offsets roughly scaled with 1/SNR im- 
plying a sky location accuracy scaling as 1/SNR^. 
Even though we did not compute proper errorboxes 
in the sky, we estimated source localization capa- 



bilities broadly consistent with theoretical expecta- 
tions derived in [2^ under similar assumptions; 

• we overestimated the frequency of sources ap- 
proaching / — lO^^Hz. This is because massive 
systems at such high frequency significantly chirp 
during the observation time (whereas chirp was not 
allowed in our template). This means that we can 
measure / and therefore estimate the chirp mass 
and, in turn, the luminosity distance of the source. 
Although this is a very appealing prospect, we esti- 
mated on average less than 1 source with a measur- 
able / in a realistic realization of the MBH binary 
population in the Universe; 

• the MS-GA performances do not seem to be af- 
fected by unequal noise levels in different MSPs. 

Our results are encouraging, however, they were still 
obtained under a number of simplifying assumptions that 
we wish to relax in our future work. Firstly, datasets 
were still evenly sampled, with no gaps; an idealised sit- 
uation that is not going to occur in reality. Secondly, we 
just took noisy datastreams and fit for the GW sources 
only, implicitly assuming perfectly known MSP param- 
eters; any realistic detection pipeline must fit for MSP 
parameters and GW signals simultaneously. Finally, we 
still did not include the pulsar terms in our injections; 
those are likely to blend together with lower frequency 
Earth terms to bias estimated source parameters and to 
(maybe) create false positives. Only by relaxing those 
assumptions we will be able to demonstrate the effective- 
ness of our MS-GA algorithm in tackling a problem with 
realistic complexity. We plan to investigate these issues 
in the next paper of the series. We will then try to apply 
our search algorithm to raw times of arrival, carrying the 
imprint of a realistic population of MBH binaries. 
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